* You may need to install the following packages: dm88_1; coefplot; center

***********************
* Spliting the sample *
***********************
use "ppc_mi.dta", replace

* split the sample into two: Sample_select and Sample_infer
splitsample, generate(sample) nsplit(2) rseed(89315)
tab sample
keep if sample==1
save "ppc_infer.dta" , replace

use "ppc_mi.dta", replace
splitsample, generate(sample) nsplit(2) rseed(89315)
tab sample
keep if sample==2
save "ppc_select.dta", replace

***********************
* Declaring variables *
***********************
global cont cps21_yob cps21_education cps21_interest_gen_1 cps21_interest_elxn_1 cps21_lr_parties_1 cps21_lr_parties_2 cps21_lr_parties_3 cps21_lr_parties_4 cps21_lr_parties_5 cps21_lr_parties_7 cps21_lr_scale_bef_1 cps21_news_cons cps21_groups_therm_1 cps21_groups_therm_2 cps21_groups_therm_3 cps21_groups_therm_4 cps21_groups_therm_6 pes21_groups1_1 pes21_groups1_2 pes21_groups1_3 pes21_groups1_4 pes21_feminine_1 pes21_masculine_1 pes21_big5_1 pes21_big5_2 pes21_big5_3 pes21_big5_4 pes21_big5_5 pes21_big5_6 pes21_big5_7 pes21_big5_8 pes21_big5_9 pes21_big5_10  cps21_children // when categories >5

global ord cps21_demsat cps21_fed_gov_sat cps21_spend_educ cps21_spend_env cps21_spend_defence cps21_spend_rec_indi cps21_covid_liberty cps21_econ_retro cps21_econ_fed_bette cps21_imm cps21_refugees cps21_govt_confusing cps21_govt_say cps21_pol_eth cps21_volunteer cps21_own_fin_retro cps21_ownfinanc_fed cps21_own_fin_future cps21_prov_gov_sat cps21_covid_sat_1 cps21_covid_sat_2 cps21_covid_sat_3 cps21_vaccine_mandat_1 cps21_vaccine_mandat_2 cps21_vaccine_mandat_3 cps21_vaccine1  pes21_dem_sat pes21_campatt pes21_keepromises pes21_losetouch pes21_hatespeech pes21_envirojob pes21_govtcare pes21_famvalues pes21_bilingualism pes21_equalrights pes21_fitin pes21_immigjobs pes21_ab_favors  pes21_ethid_1 pes21_ethid_2 pes21_ethid_3 pes21_can_id_1 pes21_can_id_2 pes21_can_id_3 pes21_can_id_4 pes21_can_id_5 pes21_can_id_6 pes21_conf_inst1_1 pes21_conf_inst1_2 pes21_conf_inst1_3 pes21_conf_inst1_4 pes21_conf_inst2_1 pes21_conf_inst2_2 pes21_conf_inst2_3 pes21_conf_inst2_4 pes21_conf_inst2_5 pes21_conf_inst2_6 pes21_discfam pes21_partic1_1 pes21_partic1_2 pes21_partic1_3 pes21_partic1_4 pes21_partic2_1 pes21_partic2_2 pes21_partic2_3 pes21_partic2_4 pes21_partic3_1 pes21_partic3_2 pes21_partic3_3 pes21_partic3_4  pes21_womenparl pes21_populism_2 pes21_populism_3 pes21_populism_4 pes21_populism_6 pes21_populism_7 pes21_populism_8 pes21_donerm pes21_donew pes21_donegl pes21_doneqc pes21_abort2 pes21_privjobs pes21_blame pes21_inequal pes21_gap pes21_hostile2 pes21_hostile4 pes21_pos_carbon pes21_pos_energy pes21_qclang pes21_qcsol pes21_newerlife pes21_health pes21_phealth pes21_mhealth pes21_lived pes21_follow_pol cps21_comfort cps21_polknow  

global bin cps21_citizenship cps21_spoil cps21_bornin_canada pes21_partymember pes21_conversion_the pes21_stdofliving pes21_trust pes21_cc1 pes21_langQC pes21_cultureQC pes21_parents_born cps21_covidrelief pes21_rural_urban

global nominal cps21_vismin_1 cps21_vismin_2 cps21_vismin_3 cps21_vismin_4 cps21_vismin_5 cps21_vismin_6 cps21_vismin_7 cps21_vismin_8 cps21_vismin_9 cps21_vismin_10 cps21_vismin_11 cps21_vismin_12 cps21_language_1 cps21_language_2 cps21_language_3 cps21_language_4 cps21_language_5 cps21_language_6 cps21_language_7 cps21_language_8 cps21_language_9 cps21_language_10 cps21_language_11 cps21_language_12 cps21_language_13 cps21_language_14 cps21_language_15 cps21_language_16 cps21_language_17  cps21_union cps21_property_1 cps21_property_2 cps21_property_3 cps21_property_4 cps21_property_5 pes21_reason_chose_1 pes21_reason_chose_2 pes21_reason_chose_3 pes21_reason_chose_4 pes21_reason_chose_5  cps21_employment1 cps21_employment2 cps21_employment3 cps21_employment4 cps21_employment5 cps21_employment6 cps21_employment7 cps21_employment8 cps21_employment9 cps21_employment10 cps21_employment11 cps21_employment12 cps21_genderid1 cps21_genderid2 cps21_genderid3 cps21_genderid4 cps21_province_1 cps21_province_2 cps21_province_3 cps21_province_4 cps21_province_5 cps21_province_6 cps21_province_7 cps21_province_8 cps21_province_9 cps21_province_10 cps21_province_11 cps21_province_12 cps21_province_13 cps21_religion1 cps21_religion2 cps21_religion3 cps21_religion4 cps21_religion5 cps21_religion6 cps21_religion7 cps21_religion8 cps21_religion9 cps21_religion10 cps21_religion11 cps21_religion12 cps21_religion13 cps21_religion14 cps21_religion15 cps21_religion16 cps21_religion17 cps21_religion18 cps21_religion19 cps21_religion20 cps21_religion21 cps21_religion22 cps21_sexuality1 cps21_sexuality2 cps21_sexuality3 cps21_sexuality4 cps21_sexuality5 cps21_sexuality6 cps21_sexuality7 cps21_marital1 cps21_marital2 cps21_marital3 cps21_marital4 cps21_marital5 cps21_marital6

* display number of variables for each type
di wordcount("cps21_yob cps21_education cps21_interest_gen_1 cps21_interest_elxn_1 cps21_lr_parties_1 cps21_lr_parties_2 cps21_lr_parties_3 cps21_lr_parties_4 cps21_lr_parties_5 cps21_lr_parties_7 cps21_lr_scale_bef_1 cps21_news_cons cps21_groups_therm_1 cps21_groups_therm_2 cps21_groups_therm_3 cps21_groups_therm_4 cps21_groups_therm_6 pes21_groups1_1 pes21_groups1_2 pes21_groups1_3 pes21_groups1_4 pes21_feminine_1 pes21_masculine_1 pes21_big5_1 pes21_big5_2 pes21_big5_3 pes21_big5_4 pes21_big5_5 pes21_big5_6 pes21_big5_7 pes21_big5_8 pes21_big5_9 pes21_big5_10  cps21_children")

di wordcount("cps21_demsat cps21_fed_gov_sat cps21_spend_educ cps21_spend_env cps21_spend_defence cps21_spend_rec_indi cps21_covid_liberty cps21_econ_retro cps21_econ_fed_bette cps21_imm cps21_refugees cps21_govt_confusing cps21_govt_say cps21_pol_eth cps21_volunteer cps21_own_fin_retro cps21_ownfinanc_fed cps21_own_fin_future cps21_prov_gov_sat cps21_covid_sat_1 cps21_covid_sat_2 cps21_covid_sat_3 cps21_vaccine_mandat_1 cps21_vaccine_mandat_2 cps21_vaccine_mandat_3 cps21_vaccine1  pes21_dem_sat pes21_campatt pes21_keepromises pes21_losetouch pes21_hatespeech pes21_envirojob pes21_govtcare pes21_famvalues pes21_bilingualism pes21_equalrights pes21_fitin pes21_immigjobs pes21_ab_favors  pes21_ethid_1 pes21_ethid_2 pes21_ethid_3 pes21_can_id_1 pes21_can_id_2 pes21_can_id_3 pes21_can_id_4 pes21_can_id_5 pes21_can_id_6 pes21_conf_inst1_1 pes21_conf_inst1_2 pes21_conf_inst1_3 pes21_conf_inst1_4 pes21_conf_inst2_1 pes21_conf_inst2_2 pes21_conf_inst2_3 pes21_conf_inst2_4 pes21_conf_inst2_5 pes21_conf_inst2_6 pes21_discfam pes21_partic1_1 pes21_partic1_2 pes21_partic1_3 pes21_partic1_4 pes21_partic2_1 pes21_partic2_2 pes21_partic2_3 pes21_partic2_4 pes21_partic3_1 pes21_partic3_2 pes21_partic3_3 pes21_partic3_4  pes21_womenparl pes21_populism_2 pes21_populism_3 pes21_populism_4 pes21_populism_6 pes21_populism_7 pes21_populism_8 pes21_donerm pes21_donew pes21_donegl pes21_doneqc pes21_abort2 pes21_privjobs pes21_blame pes21_inequal pes21_gap pes21_hostile2 pes21_hostile4 pes21_pos_carbon pes21_pos_energy pes21_qclang pes21_qcsol pes21_newerlife pes21_health pes21_phealth pes21_mhealth pes21_lived pes21_follow_pol cps21_comfort cps21_polknow")

di wordcount("cps21_citizenship cps21_spoil cps21_bornin_canada pes21_partymember pes21_conversion_the pes21_stdofliving pes21_trust pes21_cc1 pes21_langQC pes21_cultureQC pes21_parents_born cps21_covidrelief pes21_rural_urban")

di wordcount("cps21_vismin_1 cps21_vismin_2 cps21_vismin_3 cps21_vismin_4 cps21_vismin_5 cps21_vismin_6 cps21_vismin_7 cps21_vismin_8 cps21_vismin_9 cps21_vismin_10 cps21_vismin_11 cps21_vismin_12 cps21_language_1 cps21_language_2 cps21_language_3 cps21_language_4 cps21_language_5 cps21_language_6 cps21_language_7 cps21_language_8 cps21_language_9 cps21_language_10 cps21_language_11 cps21_language_12 cps21_language_13 cps21_language_14 cps21_language_15 cps21_language_16 cps21_language_17  cps21_union cps21_property_1 cps21_property_2 cps21_property_3 cps21_property_4 cps21_property_5 pes21_reason_chose_1 pes21_reason_chose_2 pes21_reason_chose_3 pes21_reason_chose_4 pes21_reason_chose_5  cps21_employment1 cps21_employment2 cps21_employment3 cps21_employment4 cps21_employment5 cps21_employment6 cps21_employment7 cps21_employment8 cps21_employment9 cps21_employment10 cps21_employment11 cps21_employment12 cps21_genderid1 cps21_genderid2 cps21_genderid3 cps21_genderid4 cps21_province_1 cps21_province_2 cps21_province_3 cps21_province_4 cps21_province_5 cps21_province_6 cps21_province_7 cps21_province_8 cps21_province_9 cps21_province_10 cps21_province_11 cps21_province_12 cps21_province_13 cps21_religion1 cps21_religion2 cps21_religion3 cps21_religion4 cps21_religion5 cps21_religion6 cps21_religion7 cps21_religion8 cps21_religion9 cps21_religion10 cps21_religion11 cps21_religion12 cps21_religion13 cps21_religion14 cps21_religion15 cps21_religion16 cps21_religion17 cps21_religion18 cps21_religion19 cps21_religion20 cps21_religion21 cps21_religion22 cps21_sexuality1 cps21_sexuality2 cps21_sexuality3 cps21_sexuality4 cps21_sexuality5 cps21_sexuality6 cps21_sexuality7 cps21_marital1 cps21_marital2 cps21_marital3 cps21_marital4 cps21_marital5 cps21_marital6")


*******************************************
* Finding the key determinants (5 chains) *
*******************************************
*** Chain 1 ***
use "ppc_select.dta", replace
keep _1_* 
foreach x of varlist _1_cps21_covid_sat_1-_1_pes21_parents_born{
	renvars `x', predrop(3)
}
lasso logit _1_ppc  c.$cont c.$ord c.$bin c.$nominal, rseed(89315)
estimates store cv
lassoknots, display(nonzero bic)

* minbic
lassoselect id=25
estimates store minbic1
lassocoef minbic1, sort(coef, standardized) nofvlabel



*** Chain 2 ***
use "ppc_select.dta", replace
keep _2_* 
foreach x of varlist _2_cps21_covid_sat_1-_2_pes21_parents_born{
	renvars `x', predrop(3)
}
lasso logit _2_ppc  c.$cont c.$ord c.$bin c.$nominal, rseed(89315)
estimates store cv
lassoknots, display(nonzero bic)

* minbic
lassoselect id=27
estimates store minbic2
lassocoef minbic2, sort(coef, standardized) nofvlabel


*** Chain 3 ***
use "ppc_select.dta", replace
keep _3_* 
foreach x of varlist _3_cps21_covid_sat_1-_3_pes21_parents_born{
	renvars `x', predrop(3)
}

lasso logit _3_ppc  c.$cont c.$ord c.$bin c.$nominal, rseed(89315)
estimates store cv
lassoknots, display(nonzero bic)

* minbic
lassoselect id=29
estimates store minbic3
lassocoef minbic3, sort(coef, standardized) nofvlabel


*** Chain 4 ***
use "ppc_select.dta", replace
keep _4_* 
foreach x of varlist _4_cps21_covid_sat_1-_4_pes21_parents_born{
	renvars `x', predrop(3)
}

lasso logit _4_ppc  c.$cont c.$ord c.$bin c.$nominal, rseed(89315)
estimates store cv
lassoknots, display(nonzero bic)

* minbic
lassoselect id=24
estimates store minbic4
lassocoef minbic4, sort(coef, standardized) nofvlabel


*** Chain 5 ***
use "ppc_select.dta", replace
keep _5_* 
foreach x of varlist _5_cps21_covid_sat_1-_5_pes21_parents_born{
	renvars `x', predrop(3)
}
lasso logit _5_ppc  c.$cont c.$ord c.$bin c.$nominal, rseed(89315)
estimates store cv
lassoknots, display(nonzero bic)

* minbic
lassoselect id=27
estimates store minbic5
lassocoef minbic5, sort(coef, standardized) nofvlabel



**************************************
* Inference of selected determinants *
**************************************
use "ppc_infer.dta", replace
set scheme white_ptol

*------------------------------*
* Standard Logistic Regression *
*------------------------------*
eststo logit_infer: mi estimate, post : logit ppc cps21_vaccine_mandat_1 cps21_vaccine_mandat_3 cps21_covid_liberty cps21_covid_sat_2  cps21_vaccine1  ///  
	pes21_conf_inst1_3 pes21_populism_2  pes21_losetouch ///
	cps21_lr_parties_2 ///
	pes21_reason_chose_1 pes21_reason_chose_5 ///
	pes21_partic3_4 pes21_partic2_1 pes21_partic1_2 ///
	pes21_pos_carbon cps21_spend_env pes21_hatespeech  pes21_gap 
estadd local cstat = "0.952", replace


* AUC calculation
preserve
mi convert flong
mi describe
local mtotal = r(M)
local cstat = 0
forvalues i = 1 / `mtotal' {
    logit ppc cps21_vaccine_mandat_1 cps21_vaccine_mandat_3 cps21_covid_liberty cps21_covid_sat_2  cps21_vaccine1  ///  
	pes21_conf_inst1_3 pes21_populism_2  pes21_losetouch ///
	cps21_lr_parties_2 pes21_reason_chose_1 pes21_reason_chose_5 ///
	pes21_partic3_4 pes21_partic2_1 pes21_partic1_2 ///
	pes21_pos_carbon cps21_spend_env pes21_hatespeech pes21_gap if _mi_m == `i'
 lroc, nog
	local cstat=`cstat'+r(area)
}
di `cstat'/`mtotal' //.95204423
restore

* coefplot
coefplot, drop(_cons)  xline(0) legend(off) ///
	coeflabels(cps21_vaccine_mandat_1="Vaccine mandate when working in a hospital" cps21_vaccine_mandat_3="Vaccine mandate when travelling by air/rail"  cps21_covid_sat_2="Dissatisfied with how provincial government handled COVID outbreak" cps21_covid_liberty="Liberty threatened by public health recommendations" cps21_vaccine1="Having not been vaccinated"  ///
	pes21_conf_inst1_3="Distrust the media"  pes21_populism_2="Compromise in politics is just selling out on one's principles"  pes21_losetouch="Those elected to Parliament soon lose touch with the people" ///
	cps21_lr_parties_2="Left-right position of CPC" ///
	pes21_reason_chose_1="* Reason for choosing the party: Party policies" pes21_reason_chose_5="* Reason for choosing the party: Other" /// 
	pes21_partic3_4="Used social media to discuss politics or political issues" pes21_partic2_1="Followed any elected officials or candidates for office on social media" pes21_partic1_2="Attended a rally or participated in a protest or demonstration" ///
	pes21_pos_carbon="The federal government should continue the carbon tax" cps21_spend_env="Spending on the environment" pes21_hatespeech ="Illegal to say hateful things publicly about racial, ethnic and religious groups"  ///	
	pes21_gap="More/less should be done to reduce the inequality gap in Canada") ///
	headings(cps21_lr_parties_2="{bf:Left/right placement}" cps21_vaccine_mandat_1="{bf:Pandemic-related variables}" ///
	pes21_conf_inst1_3="{bf:Anti-establishment attitude}" pes21_pos_carbon="{bf:Cultural attitude}"  cps21_ownfinanc_fed="{bf:Economic perception}" ///
	pes21_reason_chose_1="{bf:Reason for choosing the party}" ///
	pes21_partic3_4="{bf:Political Participation}" ///
	pes21_gap="{bf:Economic attitude}") ///
	xlabel(-1(0.5)1.5) levels(95 90) ciopts(color(black) lwidth(*1 *2)) mlcolor("0 0 128") mfcolor(white) xtitle("Standardized coefficient") ///
	text(2.0 1.2 "AUC=0.952",  size(small))
graph export "coefplot.png"	, replace 


	

* ------------------------------------------------------------------------------
* Table C.3.1
esttab  logit_infer  ///
	using "logit_infer.rtf", b(2) se(2) replace ///
	scalars ("N Observations" "cstat AUC" "rvi_avg_mi Relative increases in variance" "M_mi No. of imputations") ///
	star(+ 0.10 * 0.05 ** 0.01 *** 0.001) ///
	label wide note onecell nobaselevels  noomitted nonumbers  lines varwidth(40) ///
	coeflabels(cps21_vaccine_mandat_1 "Vaccine mandate when working in a hospital" cps21_vaccine_mandat_3 "Vaccine mandate when travelling by air/rail"  cps21_covid_sat_2 "Dissatisfied with how provincial government handled covid outbreak" cps21_covid_liberty "Liberty threatened by public health recommendations" cps21_vaccine1 "Having been vaccinated"  ///
	pes21_conf_inst1_3 "Distrust the media"  pes21_populism_2 "Compromise in politics is just selling out on one's principles"  pes21_losetouch "Those elected to Parliament soon lose touch with the people" ///
	cps21_lr_parties_2 "Left-right position of CPC" ///
	pes21_reason_chose_1 "Reason for choosing the party: Party policies" pes21_reason_chose_5 "Reason for choosing the party: Other" /// 
	pes21_partic3_4 "Used social media to discuss politics or political issues" pes21_partic2_1 "Followed any elected officials or candidates for office on social media" pes21_partic1_2 "Attended a rally or participated in a protest or demonstration" ///
	pes21_pos_carbon "The federal government should continue the carbon tax" cps21_spend_env "Spending on the environment" pes21_hatespeech  "Illegal to say hateful things publicly about racial, ethnic and religious groups" pes21_womenparl "Have more women in Parliament" ///	
	pes21_gap "More/less should be done to reduce the inequality gap in Canada") ///
	mtitles("DV: PPC voting") modelwidth(8) ///
	note ("Note: Entries are coefficients of the binomial logit model. Coefficients and standard errors are calculated according to Rubin's rules. Coefficients are standardized and standard errors are shown in parentheses.") ///
	title ("Table C.1.1 Binomial logistic regression of the selected variables") 


*------------------*
*  AUC comparisons *
*------------------*
* pandemic-related factors
preserve
mi convert flong
mi describe
local mtotal = r(M)
local cstat = 0
forvalues i = 1 / `mtotal' {
    logit ppc cps21_vaccine_mandat_3 cps21_covid_liberty cps21_covid_sat_2  ///  
	 if _mi_m == `i'
 lroc, nog
	local cstat=`cstat'+r(area)
}
di `cstat'/`mtotal' //.92262224
restore
display .92262224-.95204423

* nativist attitudes
preserve
mi convert flong
mi describe
local mtotal = r(M)
// local r2 = 0
local cstat = 0
forvalues i = 1 / `mtotal' {
    logit ppc cps21_imm pes21_fitin pes21_immigjobs ///  
	 if _mi_m == `i'
 lroc, nog
	local cstat=`cstat'+r(area)
}
di `cstat'/`mtotal' //.74769457
restore
display .74769457-.95204423


*--------------------------------*
* Rare Event Logistic Regression *
*--------------------------------*
use "ppc_infer.dta", replace
eststo relogit_infer: mi estimate, cmdok post: relogit ppc cps21_vaccine_mandat_1 cps21_vaccine_mandat_3 cps21_covid_liberty cps21_covid_sat_2  cps21_vaccine1  /// 
	pes21_conf_inst1_3 pes21_populism_2  pes21_losetouch ///
	cps21_lr_parties_2 ///
	pes21_reason_chose_1 pes21_reason_chose_5 ///
	pes21_partic3_4 pes21_partic2_1 pes21_partic1_2 ///
	pes21_pos_carbon cps21_spend_env pes21_hatespeech  pes21_gap 


coefplot, drop(_cons)  xline(0) legend(off) ///
	coeflabels(cps21_vaccine_mandat_1="Vaccine mandate when working in a hospital" cps21_vaccine_mandat_3="Vaccine mandate when travelling by air/rail"  cps21_covid_sat_2="Dissatisfied with how provincial government handled COVID outbreak" cps21_covid_liberty="Liberty threatened by public health recommendations" cps21_vaccine1="Having not been vaccinated"  ///
	pes21_conf_inst1_3="Distrust the media"  pes21_populism_2="Compromise in politics is just selling out on one's principles"  pes21_losetouch="Those elected to Parliament soon lose touch with the people" ///
	cps21_lr_parties_2="Left-right position of CPC" ///
	pes21_reason_chose_1="* Reason for choosing the party: Party policies" pes21_reason_chose_5="* Reason for choosing the party: Other" /// 
	pes21_partic3_4="Used social media to discuss politics or political issues" pes21_partic2_1="Followed any elected officials or candidates for office on social media" pes21_partic1_2="Attended a rally or participated in a protest or demonstration" ///
	pes21_pos_carbon="The federal government should continue the carbon tax" cps21_spend_env="Spending on the environment" pes21_hatespeech ="Illegal to say hateful things publicly about racial, ethnic and religious groups"  ///	
	pes21_gap="More/less should be done to reduce the inequality gap in Canada") ///
	headings(cps21_lr_parties_2="{bf:Left/right placement}" cps21_vaccine_mandat_1="{bf:Pandemic-related variables}" ///
	pes21_conf_inst1_3="{bf:Anti-establishment attitude}" pes21_pos_carbon="{bf:Cultural attitude}"  cps21_ownfinanc_fed="{bf:Economic perception}" ///
	pes21_reason_chose_1="{bf:Reason for choosing the party}" ///
	pes21_partic3_4="{bf:Political Participation}" ///
	pes21_gap="{bf:Economic attitude}") ///
	xlabel(-1(0.5)1.5) levels(95 90) ciopts(color(black) lwidth(*1 *2)) mlcolor("0 0 128") mfcolor(white) xtitle("Standardized coefficient") 
graph export "coefplot_relogit.png"	, replace 

* ------------------------------------------------------------------------------
* Table C.3.2
esttab  relogit_infer  ///
	using "relogit_infer.rtf", b(2) se(2) replace ///
	scalars ("N Observations"  "rvi_avg_mi Relative increases in variance" "M_mi No. of imputations") ///
	star(+ 0.10 * 0.05 ** 0.01 *** 0.001) ///
	label wide note onecell nobaselevels  noomitted nonumbers  lines varwidth(40) ///
	coeflabels(cps21_vaccine_mandat_1 "Vaccine mandate when working in a hospital" cps21_vaccine_mandat_3 "Vaccine mandate when travelling by air/rail"  cps21_covid_sat_2 "Dissatisfied with how provincial government handled covid outbreak" cps21_covid_liberty "Liberty threatened by public health recommendations" cps21_vaccine1 "Having been vaccinated"  ///
	pes21_conf_inst1_3 "Distrust the media"  pes21_populism_2 "Compromise in politics is just selling out on one's principles"  pes21_losetouch "Those elected to Parliament soon lose touch with the people" ///
	cps21_lr_parties_2 "Left-right position of CPC" ///
	pes21_reason_chose_1 "Reason for choosing the party: Party policies" pes21_reason_chose_5 "Reason for choosing the party: Other" /// 
	pes21_partic3_4 "Used social media to discuss politics or political issues" pes21_partic2_1 "Followed any elected officials or candidates for office on social media" pes21_partic1_2 "Attended a rally or participated in a protest or demonstration" ///
	pes21_pos_carbon "The federal government should continue the carbon tax" cps21_spend_env "Spending on the environment" pes21_hatespeech  "Illegal to say hateful things publicly about racial, ethnic and religious groups" pes21_womenparl "Have more women in Parliament" ///	
	pes21_gap "More/less should be done to reduce the inequality gap in Canada") ///
	mtitles("DV: PPC voting") modelwidth(8) ///
	note ("Note: Entries are coefficients of the binomial rare event logit model. Coefficients and standard errors are calculated according to Rubin's rules. Coefficients are standardized and standard errors are shown in parentheses.") ///
	title ("Table C.2 Binomial rare event logistic regression of the selected variables") 

		
	
********************************************
* Swapping sample_infer and sample_predict *
********************************************
*-----------------------------------------*
* Finding the key determinants (5 chains) *
*-----------------------------------------*
*** Chain 1 (swap) ***
use "ppc_infer.dta", replace
keep _1_* 
foreach x of varlist _1_cps21_covid_sat_1-_1_pes21_parents_born{
	renvars `x', predrop(3)
}
lasso logit _1_ppc  c.$cont c.$ord c.$bin c.$nominal, rseed(89315)
estimates store cv
lassoknots, display(nonzero bic)

* minbic
lassoselect id=22
estimates store minbic1
lassocoef minbic1, sort(coef, standardized) nofvlabel



*** Chain 2 (swap) ***
use "ppc_infer.dta", replace
keep _2_* 
foreach x of varlist _2_cps21_covid_sat_1-_2_pes21_parents_born{
	renvars `x', predrop(3)
}
lasso logit _2_ppc  c.$cont c.$ord c.$bin c.$nominal, rseed(89315)
estimates store cv
lassoknots, display(nonzero bic)

* minbic
lassoselect id=21
estimates store minbic2
lassocoef minbic2, sort(coef, standardized) nofvlabel


*** Chain 3 (swap) ***
use "ppc_infer.dta", replace
keep _3_* 
foreach x of varlist _3_cps21_covid_sat_1-_3_pes21_parents_born{
	renvars `x', predrop(3)
}
lasso logit _3_ppc  c.$cont c.$ord c.$bin c.$nominal, rseed(89315)
estimates store cv
lassoknots, display(nonzero bic)

* minbic
lassoselect id=28
estimates store minbic3
lassocoef minbic3, sort(coef, standardized) nofvlabel


*** Chain 4 (swap) ***
use "ppc_infer.dta", replace
keep _4_* 
foreach x of varlist _4_cps21_covid_sat_1-_4_pes21_parents_born{
	renvars `x', predrop(3)
}
lasso logit _4_ppc  c.$cont c.$ord c.$bin c.$nominal, rseed(89315)
estimates store cv
lassoknots, display(nonzero bic)

* minbic
lassoselect id=26
estimates store minbic4
lassocoef minbic4, sort(coef, standardized) nofvlabel


*** Chain 5 (swap) ***
use "ppc_infer.dta", replace
keep _5_* 
foreach x of varlist _5_cps21_covid_sat_1-_5_pes21_parents_born{
	renvars `x', predrop(3)
}
lasso logit _5_ppc  c.$cont c.$ord c.$bin c.$nominal, rseed(89315)
estimates store cv
lassoknots, display(nonzero bic)

* minbic
lassoselect id=27
estimates store minbic5
lassocoef minbic5, sort(coef, standardized) nofvlabel



*-------------------------------------------*
* Inference of selected determinants (swap) *
*-------------------------------------------*
use "ppc_select.dta", replace
set scheme white_ptol
eststo logit_infer_s: mi estimate, post: logit ppc cps21_vaccine_mandat_1 cps21_vaccine_mandat_3 cps21_covid_liberty cps21_covid_sat_1  cps21_covid_sat_3 cps21_vaccine1  /// 
	pes21_conf_inst1_3  cps21_demsat  cps21_fed_gov_sat ///
	cps21_lr_parties_2 ///
	pes21_reason_chose_1 ///
	pes21_partic1_2 pes21_partic1_4  ///
	cps21_spend_env pes21_cc1 pes21_hatespeech    //
estadd local cstat = "0.955", replace


* AUC calculation
preserve
mi convert flong
mi describe
local mtotal = r(M)
local cstat = 0
forvalues i = 1 / `mtotal' {
    logit ppc cps21_vaccine_mandat_1 cps21_vaccine_mandat_3 cps21_covid_liberty cps21_covid_sat_1  cps21_covid_sat_3 cps21_vaccine1  /// 
	pes21_conf_inst1_3  cps21_demsat  cps21_fed_gov_sat ///
	cps21_lr_parties_2 ///
	pes21_reason_chose_1 ///
	pes21_partic1_2 pes21_partic1_4  ///
	cps21_spend_env pes21_cc1 pes21_hatespeech if _mi_m == `i'
    lroc, nog
	local cstat=`cstat'+r(area)
}
di `cstat'/`mtotal' //.95529564
restore
	
* coefplot (swap)
coefplot, drop(_cons)  xline(0) legend(off) ///
	coeflabels(cps21_vaccine_mandat_1="Vaccine mandate when working in a hospital" cps21_vaccine_mandat_3="Vaccine mandate when travelling by air/rail"  cps21_covid_sat_1="Dissatisfied with how federal government handled COVID outbreak"  cps21_covid_sat_3="Dissatisfied with how local government handled COVID outbreak" cps21_covid_liberty="Liberty threatened by public health recommendations" cps21_vaccine1="Having been vaccinated"  ///
	pes21_conf_inst1_3="Distrust the media"  pes21_conf_inst1_1="Distrust the federal government"  cps21_demsat="Dissatisfied with Canadian democracy"  cps21_fed_gov_sat="Dissatisfied with the performance of the federal government"   ///
	cps21_lr_parties_2="Left-right position of CPC" ///
	pes21_reason_chose_1="* Reason for choosing the party: Party policies" pes21_reason_chose_5="Reason for choosing the party: Other" /// 
	pes21_partic1_2="Attended a rally or participated in a protest or demonstration" pes21_partic1_4="Signed a petition in person or online" ///
	pes21_pos_carbon="The federal government should continue the carbon tax" cps21_spend_env="Spending on the environment" pes21_cc1="Climate change is happening" pes21_hatespeech ="Illegal to say hateful things publicly about racial, ethnic and religious groups" pes21_abort2="Abortion should not be banned" pes21_doneqc="Less should be done for Quebec") ///
	headings(cps21_lr_parties_2="{bf:Left/right placement}" cps21_vaccine_mandat_1="{bf:Pandemic-related variables}" ///
	pes21_conf_inst1_3="{bf:Anti-establishment attitude}" cps21_spend_env="{bf:Cultural attitude}"  cps21_ownfinanc_fed="{bf:Economic perception}" ///
	pes21_reason_chose_1="{bf:Reason for choosing the party}" ///
	pes21_partic1_2="{bf:Political Participation}") ///
	xlabel(-1(0.5)1.5) levels(95 90) ciopts(color(black) lwidth(*1 *2)) mlcolor("0 0 128") mfcolor(white) xtitle("Standardized coefficient") ///
	text(2.0 1.2 "AUC=0.955",  size(small))
graph export "coefplot_swap.png", replace 



* ------------------------------------------------------------------------------
* Table C.3.2
esttab  logit_infer_s  ///
	using "logit_infer_s.rtf", b(2) se(2) replace ///
	scalars ("N Observations" "cstat AUC" "rvi_avg_mi Relative increases in variance" "M_mi No. of imputations") ///
	star(+ 0.10 * 0.05 ** 0.01 *** 0.001) ///
	label wide note onecell nobaselevels  noomitted nonumbers  lines varwidth(40) ///
	coeflabels(cps21_vaccine_mandat_1 "Vaccine mandate when working in a hospital" cps21_vaccine_mandat_3 "Vaccine mandate when travelling by air/rail"  cps21_covid_sat_1 "Dissatisfied with how federal government handled covid outbreak" cps21_covid_sat_2 "Dissatisfied with how provincial government handled covid outbreak" cps21_covid_sat_3 "Dissatisfied with how local government handled covid outbreak" cps21_covid_liberty "Liberty threatened by public health recommendations" cps21_vaccine1 "Having been vaccinated"  ///
	pes21_conf_inst1_3 "Distrust the media"  pes21_conf_inst1_1 "Distrust the federal government"  cps21_demsat "Dissatisfied with Canadian democracy"  cps21_fed_gov_sat "Dissatisfied with the performance of the federal government"   ///
	cps21_lr_parties_2 "Left-right position of CPC" ///
	pes21_reason_chose_1 "Reason for choosing the party: Party policies" pes21_reason_chose_5 "Reason for choosing the party: Other" /// 
	pes21_partic1_2 "Attended a rally or participated in a protest or demonstration" pes21_partic1_4 "Signed a petition in person or online" ///
	pes21_pos_carbon "The federal government should continue the carbon tax" cps21_spend_env "Spending on the environment" pes21_cc1 "Climate change is happening" pes21_hatespeech  "Illegal to say hateful things publicly" pes21_abort2 "Abortion should not be banned" pes21_doneqc "Less should be done for Quebec") ///
	mtitles("DV: PPC voting") modelwidth(8) ///
	note ("Note: Entries are coefficients of the binomial logit model. Coefficients and standard errors are calculated according to Rubin's rules. Coefficients are standardized and standard errors are shown in parentheses.") ///
	title ("Table C.3.2 Binomial logistic regression of the selected variables from the MI-LASSO model (Switching Sample_select and Sample_infer)") 


********************************
* Including missing indicators *
********************************
use "ppc_mi.dta", replace
* identifying which missing variables is correlated with PPC voting (regress missing indicators on PPC dummy)
foreach x of varlist cps21_covid_sat_1-pes21_parents_born {
	reg m_`x' ppc 
}

matrix m_ppc = J(274, 2, .)
local i=1
foreach x of varlist m_cps21_covid_sat_1-m_pes21_parents_born {
	reg `x' ppc 
	local coef = e(b)[1,1]
	test ppc==0
	local p = r(p)
	matrix m_ppc[`i',1] = (`coef', `p')
	local i = `i'+1
}

unab rownames: cps21_covid_sat_1-pes21_parents_born 
matrix rownames m_ppc = `rownames'
matrix list m_ppc
asdoc wmat, mat(m_ppc) replace // first column is coefficient; second column is p-value


* run lasso that includes missing indicators
use "ppc_select.dta", replace
keep _3_* m_cps21_lr_parties_5 m_cps21_lr_parties_7  m_cps21_vaccine1 m_pes21_cc1 m_pes21_groups1_2 m_cps21_econ_fed_bette m_cps21_education m_cps21_econ_retro m_pes21_groups1_3 m_pes21_masculine_1 m_pes21_pos_energy m_cps21_covid_sat_3 m_pes21_qclang m_cps21_news_cons m_pes21_groups1_4 m_pes21_qcsol m_cps21_lr_parties_1 m_cps21_lr_parties_3 m_pes21_populism_7 m_pes21_pos_carbon m_cps21_lr_parties_2 m_cps21_lr_parties_4 m_cps21_fed_gov_sat m_pes21_gap m_cps21_own_fin_future m_pes21_feminine_1 m_cps21_spend_rec_indi m_cps21_govt_say m_pes21_populism_2 m_pes21_ethid_1

foreach x of varlist _3_cps21_covid_sat_1-_3_pes21_parents_born{
	renvars `x', predrop(3)
}

* declare missing indicators
global miss m_cps21_lr_parties_5 m_cps21_lr_parties_7  m_cps21_vaccine1 m_pes21_cc1 m_pes21_groups1_2 m_cps21_econ_fed_bette m_cps21_education m_cps21_econ_retro m_pes21_groups1_3 m_pes21_masculine_1 m_pes21_pos_energy m_cps21_covid_sat_3 m_pes21_qclang m_cps21_news_cons m_pes21_groups1_4 m_pes21_qcsol m_cps21_lr_parties_1 m_cps21_lr_parties_3 m_pes21_populism_7 m_pes21_pos_carbon m_cps21_lr_parties_2 m_cps21_lr_parties_4 m_cps21_fed_gov_sat m_pes21_gap m_cps21_own_fin_future m_pes21_feminine_1 m_cps21_spend_rec_indi m_cps21_govt_say m_pes21_populism_2 m_pes21_ethid_1

* standardize the missing indicators
center $miss, addtolabel("")  standardize 
drop $miss
foreach x of varlist c_m_cps21_lr_parties_5-c_m_pes21_ethid_1{
	renvars `x', predrop(2)
}


* lasso model
lasso logit _3_ppc  c.$cont c.$ord c.$bin c.$nominal c.$miss, rseed(89315)
estimates store cv
lassoknots, display(nonzero bic)

* minbic
lassoselect id=25
estimates store minbic_miss
lassocoef minbic_miss, sort(coef, standardized) nofvlabel